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Abstract —This paper presents a fast spectral unmixing al¬ 
gorithm based on Dykstra’s alternating projection. The pro¬ 
posed algorithm formulates the fully constrained least squares 
optimization problem associated with the spectral unmixing 
task as an unconstrained regression problem followed by a 
projection onto the intersection of several closed convex sets. 
This projection is achieved by iteratively projecting onto each 
of the convex sets individually, following Dyktra’s scheme. The 
sequence thus obtained is guaranteed to converge to the sought 
projection. Thanks to the preliminary matrix decomposition and 
variable substitution, the projection is implemented intrinsically 
in a subspace, whose dimension is very often much lower than 
the number of bands. A benefit of this strategy is that the 
order of the computational complexity for each projection is 
decreased from quadratic to linear time. Numerical experiments 
considering diverse spectral unmixing scenarios provide evidence 
that the proposed algorithm competes with the state-of-the-art, 
namely when the number of endmembers is relatively small, a 
circumstance often observed in real hyperspectral applications. 

Index Terms —spectral unmixing, fully constrained least 
squares, projection onto convex sets, Dykstra’s algorithm 

1. Introduction 

S PECTRAL unmixing (SU) aims at decomposing a set 
of n multivariate measurements X = [xi,... ,Xn] into 
a collection of m elementary signatures E = [ei, • • • ,e^], 
usually referred to as endmembers, and estimating the relative 
proportions A = [ai,..., a^] of these signatures, called abun¬ 
dances. SU has been advocated as a relevant multivariate anal¬ 
ysis technique in various applicative areas, including remote 
sensing d, planetology d, microscopy d, spectroscopy 
in and gene expression analysis ID In particular, it has 
demonstrated a great interest when analyzing multi-band (e.g., 
hyperspectral) images, for instance for pixel classification 0, 
material quantification Q and subpixel detection d. 

In this context, several models have been proposed in the 
literature to properly describe the physical process underly¬ 
ing the observed measurements. Under some generally mild 
assumptions 0, these measurements are supposed to result 
from linear combinations of the elementary spectra, according 


to the popular linear mixing model (LMM) m-iH. More 
precisely, each column Xj G of the measurement matrix 
X = [xi,, Xn] can be regarded as a noisy linear combina¬ 
tion of the spectral signatures leading to the following matrix 
formulation 

X = EA + N (1) 

where 

. E G ^nxxm endmember matrix whose columns 

ei, • • • , are the signatures of the m materials, 

• A G is the abundance matrix whose jth column 

aj G contains the fractional abundances of the jth 
spectral vector Xj, 

• N G is the additive noise matrix. 

As the mixing coefficient aij represents the proportion (or 
probability of occurrence) of the the ith endmember in the 
jth measurement m, d, the abundance vectors satisfy 
the following abundance non-negativity constraint (ANC) and 
abundance sum-to-one constraint (ASC) 

Oj > 0 and = 1, Vj = 1, • • • , n (2) 

where > means element-wise greater or equal and G 
M’^xi represents a vector with all ones. Accounting for all 
the image pixels, the constraints Q can be rewritten in matrix 
form 

A>0 and l^A = 1^. (3) 

Unsupervised linear SU boils down to estimating the end- 
member matrix E and abundance matrix A from the mea¬ 
surements X following the LMM ([T]). It can be regarded as 
a special instance of (constrained) blind source separation, 
where the endmembers are the sources CD- There already 
exists a lot of algorithms for solving SU (the interested reader 
is invited to consult na-ini for comprehensive reviews 
on the SU problem and existing unmixing methods). Most 
of the unmixing techniques tackle the SU problem into two 
successive steps. Lirst, the endmember signatures are identified 
thanks to a prior knowledge regarding the scene of interest, or 
extracted from the data directly using dedicated algorithms, 
such as N-FINDR CD, vertex component analysis (VGA) 
HD, and successive volume maximization (SVMAX) ca. 
Then, in a second step, called inversion or supervised SU, the 
abundance matrix A is estimated given the previously identi¬ 
fied endmember matrix E, which is the problem addressed in 
this paper. 

Numerous inversion algorithms have been developed in 
the literature, mainly based on deterministic or statistical ap- 
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proaches. Heinz et al ini developed a fully constrained least 
squares (FCLS) algorithm by generalizing the Lawson-Hanson 
non-negativity constrained least squares (NCLS) algorithm 
na. Dobigeon et al formulated the unmixing problem into 
a Bayesian framework and proposed to draw samples from 
the posterior distribution using a Markov chain Monte Carlo 
algorithm E). This simulation-based method considers the 
ANC and ASC both strictly while the computational complex¬ 
ity is significant when compared with other optimization-based 
methods. Bioucas-Dias et al developed a sparse unmixing 
algorithm by variable splitting and augmented Lagrangian 
(SUnSAL) and its constrained version (C-SUnSAL), which 
generalizes the unmixing problem by introducing spectral 
sparsity explicitly 1^ . More recently, Chouzenoux et al ED 
proposed a primal-dual interior-point optimization algorithm 
allowing for a constrained least squares (LS) estimation ap¬ 
proach and an algorithmic structure suitable for a parallel im¬ 
plementation on modern intensive computing devices such as 
graphics processing units (GPU). Heylen et al fT2\ proposed 
a new algorithm based on the Dykstra’s algorithm ll^ for 
projections onto convex sets (POCS), with runtimes that are 
competitive compared to several other techniques. 

In this paper, we follow a Dykstra’s strategy for POCS to 
solve the unmixing problem. Using an appropriate decompo¬ 
sition of the endmember matrix and a variable substitution, 
the unmixing problem is formulated as a projection onto the 
intersection of m -1- 1 convex sets (determined by ASC and 
ANC) in a subspace, whose dimension is much lower than the 
number of bands. The intersection of m + 1 convex sets is split 
into the intersection of m convex set pairs, which guarantees 
that the abundances always live in the hyperplane governed by 
ASC to accelerate the convergence of iterative projections. In 
each projection, the subspace transformation yields linear or¬ 
der (of the number of endmembers) computational operations 
which decreases the complexity greatly when compared with 
Heylen’s method E21- 

The paper is organized as follows. In Section [III we formu¬ 
late SU as a projection problem onto the intersection of convex 
sets defined in a subspace with reduced dimensionality. We 
present the proposed strategy for splitting the intersection of 
m + 1 convex sets into the intersection of m convex set pairs. 
Then, the Dykstra’s alternating projection is used to solve this 
projection problem, where each individual projection can be 
solved analytically. The convergence and complexity analysis 
of the resulting algorithm is also studied. Section |III| applies 
the proposed algorithm to synthetic and real multi-band data. 
Conclusions and future work are summarized in Section Hvl 

H. Proposed Fast Unmixing Algorithm 

In this paper, we address the problem of supervised SU, 
which consists of solving the following optimization problem 

min||X-EA|||, 

A (4) 

subject to (s.t.) A > 0 and 

where || • \\f is the Frobenius norm. As explained in the 
introduction, this problem has been considered in many ap¬ 
plications where spectral unmixing plays a relevant role. 


It is worthy to interpret this optimization problem from a 
probabilistic point of view. The quadratic objective function 
can be easily related to the negative log-likelihood function 
associated with observations X corrupted by an additive white 
Gaussian noise. Moreover, the ANC and ASC constraints can 
be regarded as a uniform distribution for aj (Vj = 1 , • • • , n) 
on the feasible region A 


j c if aj e A 
( 0 elsewhere 


(5) 


where A = {a\a > 0, l^a = 1} and c = l/vol(^). Thus, 
minimizing @ can be interpreted as maximizing the posterior 

n 

distribution of A with the prior p(A) = Yl p{aj), where 

i=i 

we have assumed the abundance vectors are a priori 
independent. In this section, we will demonstrate that the 
optimization problem 0 can be decomposed into an uncon¬ 
strained optimization, more specifically an unconstrained least 
square (LS) problem with an explicit closed form solution, 
followed by a projection step that can be efficiently achieved 
with the Dykstra’s alternating projection algorithm. 


A. Reformulating Unmixing as a Projection Problem 

Under the assumption that E has full column ranlfl, it is 
straightforward to show that the problem @ is equivalent to 

min||Y-DA||2^ 

A (6) 

s.t. A > 0 and l^A = 

where D is any m x m square matrix such that E^E = D^D 
and 

Y = (D-^)^E^X. (7) 


Since we usually have m nx, then the formulation 
® opens the door to faster solvers. Given that E^E is 
positive definite, the equation E^E = D^D has non-singular 
solutions. In this paper, we use the Cholesky decomposition to 
find a solution of that equation. Note that we have also used 
solutions based on the eigendecomposition of E^E, leading 
to very similar results. 

Defining U = DA and the problem ® can 

be transformed as 

min||Y-U||| 

u (8) 

s.t. D-^U > 0 and b^U = 1^. 


Obviously, the optimization ([81) with respect to (w.r.t.) U 
can be implemented in parallel for each spectral vector uj, 
where U = [ui, • • • , u^] and Uj is the jth column of U. In 
another words, ® can be split into n independent problems 


mm||yj -u||i 

S.t. D“^u > 0 and b^u = 1 


(9) 


where yj is the jth column of Y (Vj = 1, • • • , n). 

Recall now that the Euclidean projection of a given vector 


^This assumption is satisfied once the endmember spectral signatures are 
linearly independent. 
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V onto a closed and convex set C is defined as El 

nc(u) = argmin (||v - u ||2 + (•c(u)) (10) 

where ic{u) denotes the indicator function 

I (X) otherwise. 

Therefore, the solution Uj of (0) is the projection 
of y-j onto the intersection of convex sets JV = 
{u G : D > 0 } (associated with the initial ANC) and 
5 = {u G : b^u = 1 } (associated with the initial ASC) 
as follows 

Uj =argmm||yj-u||| + <.A^ns(u) 

= nArns(yi) 

where Uj is the jth column of matrix U. 

Remark. It is interesting to note that Y defined by (|7]) can 
also be written as\ = DAls where Als — (E^E) E^X 
is the LS estimator associated with the unconstrained coun¬ 
terpart of Therefore, Y, U and M C^S correspond to X, 
A and A, respectively, under the linear mapping induced by 
D. 

To summarize, supervised SU can be conducted following 
Algorithm [T] by first transforming the observation matrix as 

Y = (D“^)^E^X, and then looking for the projection U 
of Y onto M r\ S. Finally, the abundance matrix is easily 
recovered through the inverse linear mapping A = D“^U. 
The projection onto J\f nS is detailed in the next paragraph. 


Algorithm 1: Fast Unmixing Algorithm 
Input: X (measurements), E (endmember matrix), M, S 

// Calculate the subspace transformation D 
from the Cholesky decomposition 

E^E = D^D 

1 D ^ Chol(E^E); 

// Compute Y 

2 Y ^ D-^E^X; 

// Project Y onto JV Cl S (Algo. O 

3 U ^ U^nsiY); 

/ / Calculate the abundance 

4 A ^ D-^U; 

Output: A (abundance matrix) 


B. Dykstra ’s Projection onto M (jS 

While the matrix Y can be computed easily and efficiently 
from d?]), its projection onto J\f (hS following (fT^ is not easy 
to perform. The difficulty mainly comes from the spectral 
correlation induced by the linear mapping D in the non¬ 
negativity constraints defining J\f, which prevents any use of 
fast algorithms similar to those introduced in l^ - lZTl dedi¬ 
cated to the projection onto the canonical simplex. However, 
as this set can be regarded as m inequalities, S CM can be 


rewritten as the intersection of m sets 

m 

Sr\N= pl^nVi 

i=l 

by splitting Af into Af = A/i fl • • • fl Afm, where A/i = 
{u G : dfu > 0 } and dj represent the ith row of D“^, 
i.e., D-i = [di,.-. . Even though projecting onto 

this m-intersection is difficult, projecting onto each convex 
set 5 n A/i (i = 1 ,..., m) is easier, as it will be shown in 
paragraph III-CI Based on this remark, we propose to perform 
the projection onto S O JV using the Dykstra’s alternating 
projection algorithm, which was first proposed in 1^ . 1^ 
and has been developed to more general optimization problems 
m, eqi. More specifically, this projection is split into m it¬ 
erative projections onto each convex set 5DA/i (i = 1 ,..., m), 
following the Dykstra’s procedure described in Algorithm [J] 


Algorithm 2: Dykstra’s Projection of Y onto S OAf 

Input: Y, D, K 

// Compute b 

1 ^ l^D-i; 

// Initial Izatlon 

2 Set ^ Y, = ... = ^ 0; 

// Main iterations 

3 for /c = 1, • • • , A do 

// Projection onto SCiAfi (Algo. [3p 

4 ^ 

5 ^ ; 

6 for i = 2, • • • , m do 

// Projection onto <SnA/i (Algo. 0 

7 uf t-n 5 nA/-,(uW 

8 ^ ^; 

9 end 

10 end 

11 

Output: U UsnAf (Y) 


The motivations for projecting onto 5 D A/i are two-fold. 
First, this projection guarantees that the vectors Uj always 
satisfy the sum-to-one constraint b^Uj = 1, which implies 
that these vectors never jump out from the hyperplane S, 
and thus accelerates the convergence significantly. Second, as 
illustrated later, incorporating the constraint b^u = 1 does 
not increase the projection computational complexity, which 
means that projecting onto S n Mi is as easy as projecting 
onto A/i (for i = 1, • • • , m). The projection onto 5 fl A/i is 
described in the next paragraph. 

C. Projection onto S CAfi 

The main step of the Dykstra’s alternating procedure (Al¬ 
gorithm |2]) consists of computing the projection U* of a given 
matrix Z onto the set 5 D A/i 

u* = nsnM (Z) 

= [n5nA/’i(zi), • • • ,n5nA/’i(Zn)]- 
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Fig. 1. Illustration of the projection of z onto the set S n A/i: the set S is 
defined by the vector c G and by the vector b orthogonal to the subspace 
<S — {c}; the vector u G <S may be written as u = Vo; + c where V spans 
the subspace *S — {c} and ol G the vector z is the orthogonal 

projection of z onto *S; the vector z is the orthogonal projection of z onto 
S n A/i, which is also the orthogonal projection of z onto the set *S n A/i. 


Let z G denote a generic column of Z. The computation 
of the projection n 5 nA/'i(z) can be achieved by solving the 
following convex constrained optimization problem: 

min||z - u ||2 


U 

s.t. dfu>0 and b^u = 1. 


(13) 


To solve the optimization ([T^ , we start by removing the 
constraint b^u = 1 by an appropriate change of variables. 
Having in mind that the set 5 = {u G : b^u = 1 } is 
a hyperplane that contains the vector c = b/||b|| 2 , then that 
constraint is equivalent to u = c +Va, where a G and 

the columns of V G span the subspace S — {c} = 

{u G : b^u = 0}, of dimension (m — 1). The matrix 

V is chosen such that V^V = Im-i^ the columns of 

V are orthonormal. Fig. [T] schematizes the mentioned entities 
jointly with z, the orthogonal projection of z onto S, and z, the 
orthogonal projection of z onto Si flA/i. The former projection 
may be written as 

z = n5(z) 

= c + P(z — c) (14) 


u = c + Va in (fTSl) . we obtain the equivalent optimization 


min ||V^ (z — c) — a ||2 s.t. (V^ d^)^ a > — (d/c) (16) 

cx 

which is a projection onto a half space whose solution is 

V^di 


a* = V-* (z - c) + Ti 


llV^di 


where 


r.=max|o,-ii^4^(V^(z-c)) 

= max{0, —sfz + fi} 




with Si = Pdi/||Pdi|| 2 , fi = —dfc/||Pdi|| 2 , and we have 
used the facts that ||V^x ||2 = ||Px ||2 and V^c = 0. 
Recalling that u = c + Va, we obtain 

Z = C + VV^(z — c) + TiSi 
V y ^ ^ 

= n5(z) -\-TiSi. 

The interpretation of dnii is clear: the orthogonal projection of 
z onto 5 n A/i is obtained by first computing z = n 5 (z), i.e.. 
the projection z onto the hyperplane S, and then computing 
z = n 5 nA/'i(z), i.e.. the projection z onto the intersection 
5 n A/i. Given that 5 fl A/i C 5, then (fTTl) is, essentially, 
a consequence of a well know result: given a convex set 
contained in some subspace, then the orthogonal projection 
of any point in the convex set can be accomplished by first 
projecting orthogonally on that subspace, and then projecting 
the result on the convex set EH Ch. 5.14], 

Finally, computing U* can be conducted in parallel for 
each column of Z leading to the following matrix update rule 
summarized in Algorithm O: 

u* =n5(z) + SiTf (18) 

with rf G given by 

rf = max{0, - sf Z} 


where 



(19) 


and the operator max has to be understood in the component¬ 
wise sense 

Note that using the Karush-Kuhn-Tucker (KKT) conditions 
to solve the problem ([T^ can also lead to this exact solution, 
as described in the Appendix. 


where P = VV^ = — bb^/||b ||2 denotes the orthogonal 

projection matrix onto S — {c}. With these objects in place, 
and given z G and u G 5, we simplify the cost function 
||z — u ||2 by introducing the projection of z onto S and by 
using the Pythagorean theorem as follows: 

||z-u||^ = ||z-z||^ + ||z-u||^ 

= ||z-z||i + ||(z-c)-Va||i 
= ||z-z||i + ||V^(z-c)-a||i (15) 

where the right hand term in (O derives directly from 
(O and from the fact that V^V = Im-i- By introducing 


D. Convergence Analysis 

The convergence of the Dykstra’s projection was first proved 
in 1^ . where it was claimed that the sequences generated 
using Dykstra’s algorithm are guaranteed to converge to the 
projection of the original point onto the intersection of the 
convex sets. Its convergence rate was explored later El, 
El- We now recall the Deutsch-Hundal theorem providing 
the convergence rate of the projection onto the intersection of 
m closed half-spaces. 

Theorem 1 (Deutsch-Hundal, 13^ : Theorem 3.8). Assuming 
that X/c is the kth projected result in Dykstra ’s algorithm and 
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Algorithm 3: Projecting Z onto S OAfi 

Input: Z, b, 

// Calculate Pd^, c and fi 

1 c^b/||b||i; 

2 Pd^ ^ d^ — cb^d^; 

3 Si ^ Pdi/||Pdi||2 ; 

4 /, ^-dfc/||Pd ,||2 ; 

// Calculate 

5 rf ^ max{0, - sf Z}; 

// Project Z onto S 

6 n 5 (Z) ^ cl^ + P(Z - cl^); 

/ / Compute the final solution U* 

7 U* ^n5(Z)+Sirf; 

Output: U* 


Xoo converged point, there exist constants 0 < c < 1 

and p > 0 such that 

||Xfe-Xoo|||</ 9 c'' (20) 

for all k. 

Theorem [T] demonstrates that Dykstra’s projection has a 
linear convergence rate 01. The convergence speed depends 
on the constant c, which depends on the number of constraints 
m and the ‘angle’ between two half-spaces O. To the best 
of our knowledge, the explicit form of c only exists for m = 2 
half-spaces and its determination for m > 2 is still an open 
problem 1351 . 

E. Complexity Analysis 

To summarize, the projection onto S Cj\f can be obtained 
by iteratively projecting onto the m sets tSflA/i (i = 1 ,..., m) 
using a Dykstra’s projection scheme as described in Algorithm 
[2j The output of this algorithm converges to the projection of 
the initial point Y onto S Off. It is interesting to note that 
the quantities denoted as in Algorithm [3] needs to be 

calculated only once since the projection of Z will be itself Z 
from the second projection n^nATs • This results from the fact 
that the projection never jumps out from the hyperplane S. 

Moreover, the most computationally expensive part of the 
proposed unmixing algorithm (Algorithm [T]) is the iterative 
procedure to project onto S Off, as described in Algorithm 
[2j For each iteration, the heaviest step is the projection onto 
the intersection 5 D A/i summarized in Algorithm [S] With 
the proposed approach, this projection only requires vector 
products and sums, with a cost of 0{nm) operations, contrary 
to the 0{nnn?) computational cost of l22l . Thus, each iteration 
of Algorithm [2] has a complexity of order 0{nnn?). 

III. Experiments using Synthetic and Real Data 

This section compares the proposed unmixing algorithm 
with several state-of-the-art unmixing algorithms, i.e., FCLS 
ifTTl . SUNSAL O, IPLS Ea and APU (H. All algorithms 
have been implemented using MATLAB R2014A on a com¬ 
puter with Intel(R) Core(TM) i7-2600 CPU@3.40GHz and 
8 GB RAM. To conduct a fair comparison, they have been 


implemented in the signal subspace without using any par¬ 
allelization. These unmixing algorithms have been compared 
using the figures of merit described in Section IIII-Al Several 
experiments have been conducted using synthetic datasets 
and are presented in Section IIII-Bl Two real hyperspectral 
(HS) datasets associated with two different applications are 
considered in Section IIII-Cl The MATLAB codes and all the 
simulation results are available on the first author’s home- 



A. Performance Measures 

In what follows. At denotes the estimation of A obtained 
at time t (in seconds) for a given algorithm. Provided that the 
endmember matrix E has full column rank, the solution of 
0 is unique and all the algorithms are expected to converge 
to this unique solution, denoted as A* = Aqo (ignoring 
numerical errors). In this work, one of the state-of-the-art 
methods is run with a large number of iterations (n = 5000 in 
our experiments) to guarantee that the optimal point A* has 
been reached. 

1) Convergence Assessment: First, different solvers de¬ 
signed to compute the solution of (H have been compared 
w.r.t. the time they require to achieve a given accuracy. Thus, 
all these algorithms have been run on the same platform and 
we have evaluated the relative error (RE) between At and A* 
as a function of the computational time defined as 


REi = 


l|A«-A1||. 

I|A*||| 


2) Quality Assessment: To analyze the quality of the un¬ 
mixing results, we have also considered the normalized mean 
square error (NMSE) 


NMSEt = 


l|At-A|||. 

I|A||| 


The smaller NMSE^, the better the quality of the unmixmg. 

II A*-AI|2 

Note that NMSEoo = iiaii 2 ^ is a characteristic of the 

j__ Il-^ll F 

objective criterion 0 and not of the algorithm. 


B. Unmixing Synthetic Data 

The synthetic data is generated using endmember spectra 
selected from the United States Geological Survey (USGS) 
digital spectral librar}0. These refiectance spectra consists 
of L = 224 spectral bands from 383nm to 2508nm. To 
mitigate the impact of the intra-endmember correlation, three 
different subsets E 3 , Eio and E 20 have been built from 
this USGS library. More specifically, Ec, is an endmember 
matrix in which the angle between any two different columns 
(endmember signatures) is larger than a (in degree). Thus, the 
smaller a, the more similar the endmembers and the higher the 
conditioning number of E. For example, E 3 contains similar 
endmembers with very small variations (including scalings) of 
the same materials and E 20 contains endmembers which are 
relatively less similar. As an illustration, a random selection 

^ http://wei.perso.enseeiht.fr/ 

^ http:// speclab. cr. usgs. gov/ spectral. Iib06/ 
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of several endmembers from E 3 and E 20 have been depicted 
in Fig. [21 The abundances have been generated uniformly in 
the simplex A defined by the ANC and ASC constraints. 

Unless indicated, the performance of these algorithms has 
been evaluated on a synthetic image of size 100 x 100 whose 
signal to noise ratio (SNR) has been fixed to SNR=30dB and 
the number of considered endmembers is m = 5. 




Fig. 2. Five endmember signatures randomly selected from E 3 (left) and 
E 20 (right). 


Section [riI-B2l tSNR= 30dB and n = 100^). The endmember 
signatures have been selected from Eio (similar results have 
been observed when using E 3 and E 20 ). All the algorithms 
have been stopped once At reaches the same convergence 
criterion REt < — lOOdB. The proposed SUDAP has been 
compared with the two most competitive algorithms SUNSAL 
and APU. The final REs and the corresponding computational 
times versus m have been reported in Fig. [4l including error 
bars to monitor the stability of the algorithms (these results 
have been computed from 30 Monte Carlo runs). 




1) Initialization: The proposed SUDAP, APU and FCLS 
algorithms do not require any initialization contrary to SUN- 
SAL and IPLS. As suggested by the authors of these two 
methods, SUNSAL has been initialized with the unconstrained 
LS estimator of the abundances whereas IPLS has been 
initialized with the zero matrix. Note that our simulations have 
shown that both SUNSAL and IPLS are not sensitive to these 
initializations. 

2) Performance vs. Time: The NMSE and RE for these 
five different algorithms are displayed in Eig. [3] as a function 
of the execution time. These results have been obtained by 
averaging the outputs of 30 Monte Carlo runs. More pre¬ 
cisely, 10 randomly selected matrices for each set E 3 , Eio 
and E 20 are used to consider the different intra-endmember 
correlations. All the algorithms converge to the same solution 
as expected. However, as demonstrated in these two figures, 
SUNSAL, APU and the proposed SUDAP are much faster 
than ECLS and IPLS. Prom the zoomed version in Pig. (S) we 
can observe that in the first iterations SUDAP converges faster 
than APU and SUNSAL. More specifically, for instance, if the 
respective algorithms are stopped once RE^ < — 80dB (around 
t = 50ms), SUDAP performs faster than SUNSAL and APU 
and with a lower NMSE^. 




Fig. 3. NMSE (left) and RE (right) V 5 . time (zoomed version in top right). 

3) Time vs. the Number of Endmembers: In this test, the 
number of endmembers m varies from 3 to 23 while the 
other parameters have been fixed to the same values as in 


Fig. 4. RE (left) and time (right) v^. number of endmembers for SUNSAL, 
APU and SUDAP (REt < -lOOdB). 

Pig. in (left) shows that all the algorithms have converged to 
a point satisfying REt < — lOOdB and that SUDAP and APU 
are slightly better than SUNSAL. However, SUNSAL provides 
a smaller estimation variance leading to a more stable estima¬ 
tor. Pig. [4| (right) shows that the execution time of the three 
methods is an increasing function of the number of endmem¬ 
bers m, as expected. However, there are significant differences 
between the respective rates of increase. The execution times 
of APU and SUDAP are cubic and quadratic functions of m 
whereas SUNSAL benefits from a milder increasing rate. More 
precisely, SUDAP is faster than SUNSAL when the number 
of endmembers is small, e.g., smaller than 19 (this value may 
change depending on the SNR value, the conditioning number 
of E, the abundance statistics, etc.). Conversely, SUNSAL is 
faster than SUDAP for m > 19. SUNSAL is more efficient 
than APU for m > 15 and SUDAP is always faster than 
APU. The error bars confirm that SUNSAL offers more stable 
results than SUDAP and APU. Therefore, it can be concluded 
that the proposed SUDAP is more promising to unmix a 
multi-band image containing a reasonable number of materials, 
while SUNSAL is more efficient when considering a scenario 
containing a lot of materials. 

4) Time vs. Number of Pixels: In this test, the performance 
of the algorithms has been evaluated for a varying number of 
pixels n from 100 ^ to 400^ (the other parameters have been 
fixed the same values as in Section IIII-B21) . The endmember 
signatures have been selected from Eio (similar results have 
been observed when using E 3 and E 20 ) and the stopping 
rule has been chosen as RE^ < — lOOdB. All results have 
been averaged from 30 Monte Carlo runs. The final REs 
and the corresponding computational times are shown in Pig. 
[ 5 ] The computational time of the three algorithms increases 
approximately linearly w.r.t. the number of image pixels and 
SUDAP provides the faster solution, regardless the number of 
pixels. 













































































7 




Fig. 5. RE (left) and time (right) vs. number of pixels for SUNSAL, APU 
and SUDAP (REt < -lOOdB). 

5) Time vs. SNR: In this experiment, the SNR of the HS 
image varies from OdB to 50dB while the other parameters 
are the same as in Section IIII-B2I The stopping rule is the 
one of Section UlI-BSI The results are displayed in Fig. [6] and 
indicate that SUNSAL is more efficient than APU and SUDAP 
(i.e., uses less time) for low SNR scenarios. More specifically, 
to achieve REt < — lOOdB, SUNSAL provides more efficient 
unmixing when the SNR is lower than 5dB while SUDAP is 
faster than SUNSAL when the SNR is higher than 5dB. 




Eig. 6. RE (left) and time (right) V5. SNR for SUNSAL, APU and SUDAP 
(REt < -lOOdB). 


C. Real Data 

This section compares the performance of the proposed 
SUDAP algorithm with that of SUNSAL and APU using two 
real datasets associated with two different applications, i.e., 
spectroscopy and hyperspectral imaging. 

1) EELS Dataset: In this experiment, a spectral image 
acquired by electron energy-loss spectroscopy (EELS) is con¬ 
sidered. The analyzed dataset is a 64 x 64 pixel spectrum- 
image acquired in nx = 1340 energy channels over a region 
composed of several nanocages in a boron-nitride nanotubes 
(BNNT) sample (S). A false color image of the EELS data 
(with an arbitrary selection of three channels as RGB bands) 
is displayed in Eig. [ 7 ] (left). Eollowing 0 , the number of end- 
members has been set to m = 6. The endmember signatures 
have been extracted from the dataset using VGA Ha and are 
depicted in Eig. [ 7 ] (right). The abundance maps estimated by 
the considered unmixing algorithms are shown in Eig. 0 for a 
stopping rule defined as RE^ < lOOdB. 

There is no visual difference between the abundance maps 
provided by SUNSAL, APU and the proposed SUDAP. Since 
there is no available ground-truth for the abundances, the 
objective criterion Jt = ||X —EAt|||^ minimized by the algo¬ 
rithms has been evaluated instead of NMSE^. The variations of 



Fig. 7. EELS dataset: HS image (left) and extracted endmember signatures 
(right). 


the objective function and the corresponding REs are displayed 
in Eig. as a function of the computational time. 




Fig. 9. Objective (left) and RE (right) vs. time for SUNSAL, APU and 
SUDAP (EELS data). 

Both figures show that the proposed SUDAP performs faster 
than APU and SUNSAL as long as the stopping rule has been 
fixed as REt < — 60dB. Eor lower REt, SUDAP becomes 
less efficient than SUNSAL. To explore the convergence more 
explicitly, the number of spectral vectors that do not satisfy 
the convergence criterion, i.e., for which RE> — lOOdB, 
has been determined and is depicted in Eig. [TOl It is clear 
that most of the spectral vectors (around 3600 out of 4096 
pixels) converged quickly, e.g., in less than 0.02 seconds. The 
remaining measurements (around 500 pixels) require longer 
time to converge, which leads to the slow convergence as 
observed in Eig. 0 The slow convergence of the projection 
methods for these pixels may result from an inappropriate 
observational model due to, e.g., endmember variability ll3^ 
or nonlinearity effects im. On the contrary, SUNSAL is more 
robust to these discrepancies and converges faster for these 
pixels. This corresponds to the results shown in Eig. 0 



Fig. 10. Number of pixels that do not satisfy the stopping rule vs. time for 
SUNSAL, APU and SUDAP (EELS data). 


































































Fig. 8. EELS dataset: abundance maps estimated by SUNSAL (top), APU (middle) and SUDAP (bottom). 


2) Cuprite Dataset: This section investigates the perfor¬ 
mance of the proposed SUDAP algorithm when unmixing 
a real HS image. This image, which has received a lot 
of interest in the remote sensing and geoscience literature, 
was acquired over Cuprite field by the JPL/NASA airborne 
visible/infrared imaging spectrometer (AVIRIS) 153. Cuprite 
scene is a mining area in southern Nevada composed of 
several minerals and some vegetation, located approximately 
200km northwest of Las Vegas. The image considered in this 
experiment consists of 250 x 190 pixels of nx = 189 spectral 
bands obtained after removing the water vapor absorption 
bands. A composite color image of the scene of interest is 
shown in Fig. [TT] (left). As in Section lni-C21 the endmember 
matrix E has been learnt from the HS data using VC A. 
According to ca, the number of endmembers has been set to 
m = 14. The estimated endmember signatures are displayed 
in Fig. [TT] (right) ant the first five corresponding abundance 
maps recovered by SUNSAL, APU and SUDAP are shown in 
Fig. Visually, all three methods provide similar abundance 
mapo 

From Fig. [TT] (right), the signatures appear to be highly 
correlated, which makes the unmixing quite challenging. This 
can be confirmed by computing the smallest angle between any 
couple of endmembers, which is equal to a = 2.46 (in degree). 
This makes the projection-based methods, including SUDAP 
and APU, less efficient since alternating projections are widely 
known for their slower convergence when the convex sets 
exhibit small angles, which is consistent with the convergence 
analysis in Section III-Dl Fig. [TSj which depicts the objective 
function and the RE w.r.t. the computational times corroborates 

"^Similar results were also observed for abundance maps of the other 
endmembers. They are not shown here for brevity and are available in a 
separate technical report (38). 



Fig. 11. Cuprite dataset: HS image (left) and extracted endmember signatures 
(right). 


this point. Indeed, SUDAP performs faster than SUNSAL and 
APU if the algorithms are stopped before RE< —30dB. For 
lower REt, SUNSAL surpasses SUDAP. 



Time(s) Time(s) 


Fig. 13. Objective function (left) and RE (right) vs. time for SUNSAL, APU 
and SUDAP (Cuprite data). 


IV. Conclusion 

This paper proposed a fast unmixing method based on 
an alternating projection strategy. Formulating the spectral 
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Fig. 12. Cuprite dataset: abundance maps estimated by SUNSAL (top), APU (middle) and SUDAP (bottom). 


unmixing problem as a projection onto the intersection of 
convex sets allowed Dykstra’s algorithm to be used to compute 
the solution of this unmixing problem. The projection was 
implemented intrinsically in a subspace, making the proposed 
algorithm computationally efficient. In particular, the proposed 
unmixing algorithm showed similar performance comparing to 
state-of-the-art methods, with significantly reduced execution 
time, especially when the number of endmembers is small or 
moderate, which is often the case when analyzing conventional 
multi-band images. Future work includes the generalization of 
the proposed algorithm to cases where the endmember matrix 
is rank deficient or ill-conditioned. 


Appendix 

Solving ([T3]) with KKT conditions 

Following the KKT conditions, the problem ([13 can be 
reformulated as finding u* satisfying the following conditions 

u* — z + /ib — Adi 
dfu* 
b^u* 

A 

Adfu* 

Direct computations lead to 

u* = z — z + Az (22) 


> 

> 

> 


0 

0 

1 

0 

0 

0 . 


( 21 ) 


where 


z 

= c(b^z - 

1) 

c 

= b/||b||2 


Az 

= nsi 


Ti 

= max{0, - 

-df (z - 

Si 

= Pdi/||Pdi||2 

P 

II 

HH 

1 

cr 

cr 

7l|b||^ 


llPddb} 


( 23 ) 
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Computing the projection of zj for j = 1, • • • ,n can be 
conducted in parallel, leading to the following matrix update 
rule 

U* =n5nA/-,(Z) 

= Z - Z + SiTf (24) 

= n5(z) + Sixf 

where 

Z = c (b^Z - ll) 

rf = max{ 0 , -df - zj /||Pdi|| 2 }. 

As a conclusion, the updating rules (1^ and ([T^ only differ 
by the way the projection Ilsi'Zi) onto S has been computed. 
However, it is easy to show that fls{'^) = Z — Z used in (1^ 
is fully equivalent to n 5 (Z) = cl^ + P(Z — cl^) required 
in ([TS]) . 

Remark. It is worthy to provide an alternative geometric 
interpretation of the KKT-based solution (1^ . First, z — z 
is the projection of z onto the affine set S. Second, if the 
projection is inside the set Mi, which means df (z — z) > 0, 
then the projection of z onto the intersection S F Mi is 
z — z. If the projection is outside of the set Mi, implying that 
df (z — z) < 0 , a move Az inside the affine set S should 
be added to z — z to reach the set Mi. This move Az should 
ensure three constraints: I) Az keeps the point z — z + Az 
inside the affine set S, 2) z — z + Az is on the boundary of the 
set Mi, and 3) the Euclidean norm of Az is minimal. The first 
constraint, which can be formulated as b^Az = 0, is ensured 
by imposing a move of the form Az = Pw where P = VV^ 
is the projector onto the subspace So orthogonal to b. The 
second constraint is fulfilled when df (z — z + Az) = 0, 
leading to df Pw = —5i, where Si = df (z — z). Thus, due 
to the third constraint, w can be defined as 

w = argmin ||Pv ||2 s.t. dfPw = —Si. (25) 

V 

Using the fact that P is an idempotent matrix, i.e., P^ = P, 
the constrained optimization problem can be solved analyti¬ 
cally with the method of Lagrange multipliers, leading to 

w =-(5j(dfPdj) ^dj (26) 

and Az = Pw = —^^(dfPdi) ^Pd^. This final result is 
consistent with the move defined in ( 1 ^ and ( 1 ^ by setting 
Ti = max{Q, — ||pd^|| 2 } ~ Pdi/||Pd^|| 2 . Recall that 

||Pdi ||2 = (dfPdij since P^P = P. 
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